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NUMERICAL RESULTS FOR THE DIFFRACTION OF A NORMAL SHOCK WAVE 
BY A SPHERE AND FOR THE SUBSEQUENT TRANSIENT FLOW 

By Richard W. Barnwell 
Langley Research Center 

SUMMARY 

The finite -difference method which Peter D. Lax developed for treating unsteady 
inviscid flow fields is used to study the transient flow in the shock layer of a sphere that 
has been struck by a normal shock wave. Transient flow of this sort is encountered when 
a shock tube is used as a supersonic wind tunnel. Time histories of the shock detach- 
ment distance and the stagnation -point pressure and tangential velocity gradient are pre- 
sented for ranges of the incident-shock Mach number and the perfect-gas specific-heat 
ratio. These results show that the stagnation- point pressure approaches the steady value 
much more rapidly than the shock detachment distance. In general, the stagnation-point 
tangential velocity gradient approaches the steady value less rapidly than stagnation-point 
pressure but more rapidly than the shock detachment distance. As the specific-heat 
ratio is decreased and the incident-shock Mach number is increased, the variation of the 
velocity gradient with respect to the shock detachment distance becomes more nearly 
linear. 


INTRODUCTION 

Supersonic flow about a model mounted in a shock tube is initiated by causing a 
normal shock wave to travel down the tube past the model. When the shock strikes the 
model it is diffracted, and a shock layer is formed ahead of the model. This shock layer 
adjusts with time until steady flow is established. 

Several investigators have obtained approximate analytical expressions for the 
time history of the shock detachment distance for a blunt model which has been washed 
by a normal shock by assuming that the model is placed in a still gas and accelerated 
instantaneously to supersonic speed. Cabannes (ref. 1) has developed a second-order 
Taylor series expansion for smooth plane and axisymmetric blunt bodies. Bausset 
(ref. 2) has obtained the third- and fourth-order terms for this expansion and has 
included three-dimensional effects. Miles, Mire Is, and Wang (ref. 3) have used a mass- 
balance treatment for the special case of the flat-nosed cylinder. 



Experimental results have been obtained for the transient behavior in the flow 
fields of various blunt bodies which have been washed by normal shocks. Results have 
been presented by Davies (ref. 4) for flat-nosed cylinders with rounded corners. In 
addition, he has derived an expression for the time variation of the shock-layer thickness 
by assuming that the bow-shock deceleration is constant. Offenhartz and Weisblatt 
(ref. 5) have obtained data for circular cylinders alined with the flow with several nose 
bluntnesses. Some results for cylinders with axes perpendicular to the flow and for 
spheres have been published in references 6, 7, and 8. 

Numerical results for diffraction of normal shock waves by flat-nosed cylinders 
coaxial with the flow have been obtained by Butler (ref. 9). 

In this report numerical results are presented for the time histories of the shock 
detachment distance and the stagnation-point pressure and tangential velocity gradient 
for a sphere which has been washed by a normal shock wave. The flow field is assumed 
to be inviscid, and a calorically perfect gas model is employed. The cases which are 
treated cover a wide range of incident -shock Mach number and specific-heat ratio. A 
finite -difference method based on a technique of first-order accuracy which Lax (ref. 10) 
developed for treating one -dimensional inviscid unsteady flow fields containing shock 
waves is used to calculate the results. The Lax technique is used to make blunt-body 
flow-field calculations in references 11 and 12. Other applications of the Lax technique 
are given in references 13 and 14. 


SYMBOLS 

Ai,Bi, Ci,Di,Ei quantities defined by equations (10) 
a speed of sound 

c speed of reflected shock wave for one -dimensional shock reflection problem 

K constant defined by equations (11) 

Kl constant defined by equations (14) 

Ms incident-shock-wave Mach number 

M Mach number 

p pressure 
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r radial distance from center of sphere 

t time 

u radial velocity component 

V total magnitude of velocity 

v tangential velocity component 

a constant defined by equation (A8) 

y perfect-gas specific-heat ratio 

At,Ar,A0 mesh spacings 
6 shock detachment distance 

6 angular coordinate measured from axis of symmetry 

v constant defined by equations (9) 

p mass density 

t nondimensional time defined by equation (17) 

Subscripts: 

b at the body surface 

l spatial position, r = r^ + (l - l)Ar 

m spatial position, 6 - (m - 1)A0 

max maximum 

rs quantity evaluated behind the reflected shock at the point where a normal 

shock wave is diffracted by a surface 
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S steady flow 

s at the shock 

stag at the stagnation point 

t = 0 immediately after shock impingement 

°° free stream 

Superscript: 

k number of time steps, t = k At 

A prime denotes differentiation with respect to 6. 

A bar denotes a nondimensional quantity as given by equations (A4). 

PROBLEM DESCRIPTION 

The transient behavior which occurs when a normal shock wave moves into a still 
gas and strikes a sphere is to be investigated. The flow field is considered to be invis- 
cid, and a calorically perfect gas model is assumed. Phenomena which are treated 
include the diffraction of the incident shock wave by the sphere and the unsteady flow in 
the shock layer ahead of the body after shock impingement and before the establishment 
of steady flow. 

Numerical solutions are obtained to the initial-value problem for the partial dif- 
ferential equations which govern unsteady inviscid gas flow and the boundary conditions 
for flow over a sphere. The flow properties are determined at successive instants in 
time at points fixed in the flow field. 

Governing Equations 

The natural coordinate system for treating flow about a sphere is a spherical polar 
system with its origin at the center of the sphere. Since the flow is symmetric about the 
normal to the incident shock wave which passes through the center of the sphere, the axis 
of the coordinate system is chosen to lie along this line. The flow properties are func- 
tions of the distance from the center of the sphere r, the polar angle 6, and the time t. 

The usual conservation form of the governing equations for unsteady axisymmetric 
flow in terms of spherical polar coordinates is given in reference 11. As is pointed out 
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in that reference, this form of the governing equations cannot be used to evaluate the flow 
properties along the axis where 0 and sin 6 are equal to zero. A useful conserva- 
tion form of the equations for application along the axis is obtained by dividing the equa- 
tions in reference 11 by r sin 9, the perpendicular distance from the axis, and evaluating 
the indeterminate terms which result with the aid of l'Hospital’s rule. 

When the equations of reference 11 are divided by r sin 9, the following equations 
result: 

Continuity: 


9(pr) f 9 (pur) ^ 8 (py) 
8t 8r 8 9 


+ pu + pv cot 


0 = 0 


Radial momentum: 


8 (pur) 
8t 


( „ 2 s ) 

8 (puv) 

/ o o\ 

(j? + pu“)r 

80 {_ 

p - pfu^ - V*) 


+ puv cot 0 = 0 


( 1 ) 


( 2 ) 


Tangential momentum: 

8 ^ r '' + -^(puvr) + -^(p + pv2) + 2puv + pv 2 cot 0 = 0 (3) 

Energy: 



In this report these equations are used for all values of 0 except 0 = 0, where 
the terms containing the product v cot 0 are indeterminate. As previously mentioned, 
l'Hospital's rule is used to evaluate these indeterminate forms. The governing equations 
for use when 0 = 0 are written as: 

Continuity: 


8(p^ + 8(purl + 2 9(pvl + pu = 0 (5) 

8t 8r 80 
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Radial momentum: 


-^(Pur) + ^|(p + pu2)rj + 2 -|^(puv) 



Energy: 


_a_ 

at 




p(u 2 + y2) 





P + 1 P(u 2 + 



+ 2 


_0_ 

00 



P + 1 p(u 2 + V 2 ) 


§ +u C- 



= 0 


( 6 ) 


(7) 


The tangential momentum equation is not needed when 0 = 0 because the tangential 
velocity component is known to have a value of zero on the axis. 


Boundary Conditions 

The computational region which is used in this report is shown in figure 1. It is 
bounded by two circular arcs and two radial lines. Along the arc r = r^, which 



Figure 1.- Computational region. 
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represents the surface of the sphere, the condition u = 0 is applied since the flow 
direction must be tangent to the surface at the body. Along the axis, where 0 = 0, the 
symmetry condition v = 0 is applied. The arc r = r max is positioned so that it lies 
outside the steady-shock layer. The line 0 = 0max is positioned so that the flow at 
points in its vicinity is supersonic at all times after the incident shock wave has passed 
out of the computational region. 

COMPUTATIONAL METHOD 


Finite -Difference Equations 

The computational region shown in figure 1 is divided with a polar grid system, and 
the flow properties are computed at the grid intersections at successive instants in time. 
The mesh spacings Ar, A 0, and At are constants such that the ratios At/Ar and 
At/A0 do not exceed the limits specified by the Courant-Friedrichs-Lewy stability cri- 
teria, which are discussed later. The indices for the time t and the coordinates r 
and 0 are k, Z,and m, respectively. 

The partial differential equations (1) to (7) are represented by the index equations 


0Aj dB i 3Cs 

— — + - + (1 + v )— i + Di + (1 - y)Ei cot 0 = 0 


at ar 


30 


(i = 1,2,3, 4) 


where 


and where 


v= 1 
v = 0 


(0 = 0 ) 
(0 > 0 ) 


> 

I — 1 

II 

*0 

A2 = pur 

A3 = pvr 

Bj = pur 

B2 = (p + pu 2 )r 

B3 = puvr 

O 

h- 1 

II 

$ 

C2 = puv 

C3 = p + pv 

Dj = pu 

D2 = -jp - P(u 2 - v 2 )J 

D3 = 2 puv 

Ej = pv 

E2 = puv 

E3 = pv 2 


A4 = | p + \ p(u2 + v2 ^] ] 


E4 = 


= v — p + -i p(u 2 + v 2 )j 


( 8 ) 


(9) 


> ( 10 ) 
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The Lax technique can be extended to the problem which is being treated in this 
report by replacing the partial derivatives (9Ai/9t)^ m at the point (k,l,m) with the for- 
ward difference expressions 


where 


< A C 


1+1 ,m 


+ (Ai) 


1-1, m 


At 




K l (Ai) l,m + l + (Ai) 


•{-$ 


K = 1 

K = sec A$ 


(i = 1,4) 1 
(i = 2,3) | 


and by using the central difference expressions 


and 


(Bi) 


k 

1+1 ,m 


- (Bi) 


k 

1-1, m 


2 Ar 


( 11 ) 


(Ci) 


l,m+l 


(Ci) 


l,m-l 


2 A6 


k k 

to represent the partial derivatives (9Bj/9r) and (9Ci/9 0) , respectively, at this 

C , m. if jHl 

point. The finite -difference expressions for the time derivatives in the momentum equa- 
tions (i = 2,3) are modified so that they will vanish for uniform flow parallel to the axis. 
For uniform flow, A2 and A3 are not constant but vary as cos 6 and sin 6, respec- 
tively, and 


(A i>f,m+l + <Al) ?,m-l = 2 < Ai >*m cos iS (i = 2 > 3) 

If these sums were not modified, errors of order (A0)^ would be introduced into the 

k k 

solution for the uniform flow with each computation. The terms (Di) and (Ei) 

L,m £,m 

at the point (k,l,m) are replaced by the expressions 
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and 




(Ei)^ i + (Ei), k , + (Ei) 7 k , + (Ei) k , 

Z+ l,m v L 'l- l,m i, m+1 1; Z, m-1 


respectively. 


The finite -difference expressions described above are used to obtain finite- 
difference representations for the partial differential equations (8) at the point (k,Z,m). 
They are written as 


(Ai) 


k+1 

Z,m 


4 ^ Ai ^+l,m + + f Ai ^,m+l 



N k “I 1 At r„ .k 
- (Bi) „ + o~T2 (^i) i 



(1 + v) - i At 



k 

Z+l,m 


+ (Di) 


k 

Z-l,m 


+ 



m+1 



+ (Ei) 


k 

Z,m+1 


+ " v)coi 9 (i = 1,2, 3, 4) (12) 

All the quantities on the right-hand side of equations (12) have the time index k and can 
be evaluated with the aid of equations (10) since the flow properties are known at time k. 

Only the four points (k, Z+l,m), (k, Z-l,m), (k, Z, m+1), and (k, Z,m-1) are used to 
evaluate the right-hand side of equations (12); the point (k,Z,m) is not needed. Therefore, 
a staggered grid is used in this report, as shown in figure 2. 

When the terms in equations (12) are replaced by Taylor series expansions about 
the point (k,Z,m), the following equations result: 


/0Aj\ k / 9Bj\ k (d CA 1 

hr + tar + (1 + 4r? 

\ 8t 'Z,m 

k 




a 2 Ai 


+ m + (1 " m cot 9 

l, m l ’ m l ’ m 


(Ar) 2 ^ 2 A i 


at * 1 




4 At 




\ , i (Ae) 2 fe 2A i\ k 

I 4 At 
Z,m 


5 « 2 4, m 


i^(Ai), jm K lt 0(i 2 ) (13) 


(i = 1,2, 3, 4) 
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where 







Equations (13) differ from equations (8) because of the presence of the terms on the right- 
hand side. Since the second and third terms are similar to some of the viscous terms in 
the viscous momentum equations, they are designated as "artificial viscosity" terms. 

The presence of these terms causes the inviscid flow properties to be continuous across 
shock waves. As a result, the shocks are smeared over several mesh spaces, and com- 
putations can be performed in their vicinity in the same manner as elsewhere. It can be 
shown that for uniform flow all the terms of order A on the right-hand side of equa- 
tions (13) either vanish or cancel. The cancellation of terms in the momentum equations 
(i = 2,3) occurs because of the presence of the weighting function K in the finite- 
difference equations (12). 



1 = 5 4 3 2 1 

Figure 2.- Sample grid. 
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In this report equations (12) are used to determine the flow properties at points 
with an odd time index, and an alternate set of difference equations is used to obtain the 
flow properties at the points which have an even time index. This alternate set of equa- 
tions does not contain implicit "artificial viscosity" terms. It has been found that the 
solutions which this alternating procedure yields are stable, and that the results are 
more accurate than those obtained when equations (12) are used exclusively. 

The alternate set of finite -difference equations is obtained by replacing the quanti- 
k k k 

ties (3Ai/8t) , (Di) , and (Ei) at the point (k,Z,m) by the finite -difference 

Z,m l, m l,m 

expressions 


Km - I^L,m + (A i>f-l,m) + | < Ai) f + 2,m~ 2( H,m + (Ai) z-2,m 


At 


i[<i, 


m + 


and 


< E ‘>W,m + (Ei > 




respectively, and using the same expressions for the partial derivatives (9Bi/3r)_ 

c j m 

and (9Ci/30) as are used to obtain equations (12). When these expressions are sub- 
Z,m 

stituted into equations (8), the following finite -difference equations result: 


+ • 2( < + K», 

- ®»M,m] " < C <m + ]- I 4 *K + l, m + ®»M, m ] 


m 




(Ei) f+l,m + (Ei) Z- 


1,J 


cot e 


(i = 1,2, 3, 4) 


(15) 


Note that equations (15) are written for a staggered grid. 
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Results obtained by using equations (12) exclusively and by using equations (12) and 
(15) alternately are compared subsequently. 

Boundary Points 

In order to compute the flow properties at a point on the axis {9 = 0) with either 
equations (12) or (15), the flow properties for the previous time are required at points on 
either side of the axis. One of these points lies outside the region of computation. The 
flow properties at this outside point are evaluated by reflecting the properties at the cor- 
responding point within the region of computation. The tangential velocity component v 
is an odd function of 9, and the other properties are even. 

Equations (12) and (15) are used to determine the flow properties at points on the 
surface of the sphere as well as in the flow field. When these equations are used at the 
surface of the sphere or when equations (15) are used on the arc next to the surface, 
values are required for the flow properties at fictitious points which are located within 
the sphere. These values are obtained by reflecting the values at the points which lie an 
equal distance outside the surface along the same radial line. The radial velocity com- 
ponent u is treated as an odd function, and the other flow properties are considered to 
be even functions. One consequence of this procedure is that the finite -difference expres- 
sions which represent the normal derivatives (9Bi/9r)„_ T . (i = 1,3,4) at the body vanish. 

i-ib 

This is not physically correct, but it is tolerated since it does not cause appreciable error. 
The radial momentum equation (i = 2) is not used at the surface since the radial com- 
ponent u is known to be zero when r = r^. 

The flow properties at points along the arc r = r max are those of either the undis- 
turbed gas or the uniform flow behind the incident shock wave. During the early part of 
the computation when the incident shock wave intersects this arc, the location of the inter- 
section is determined from the incident-shock speed and the field of computation, and the 
flow properties at points on the arc are assigned accordingly. Uniform flow properties 
are maintained along the arc when the incident shock no longer intersects it. 

The incident shock wave is advanced along the line 6 = 0 max in a manner similar 
to that which is used for the arc r = r max . When the incident shock wave has passed out 
of the computational region, the flow properties along the line 6 = 0 max are determined 
by extrapolating the quantities Cj linearly. The line 9 = 6 max is located in a region 
where the flow is supersonic at all times after the incident shock wave has passed, because 
the extrapolation procedure is not stable if the flow is locally subsonic. 
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Stability Criteria 


The Courant-Friedrichs-Lewy stability conditions (ref. 15) for the partial differ- 
ential equations which govern inviscid, compressible fluid flow place an upper bound on 
the mesh spacing At for particular values of the other mesh spacings. These condi- 
tions for a spherical polar coordinate system are 


At < 1 

Ar ” |u 1+ a 


At_ < r 
A6» | v | + a 


(16) 


If these conditions are not met the solution at time k + 1 will not be influenced properly 
by the solution at previous times and unstable computations can result. The present 
research indicates that these stability conditions, with the limitation mentioned below, 
are sufficient. 


Since the mesh spacings are constant for a given computation, the smallest values 
of the terms on the right-hand side of the inequalities (16) must be used. According to 
exact calculations, the smallest value for the right-hand side of the inequality for At/Ar 
for the cases treated is (Voo + a 00 ) - '*'. However, it was found that the use of this value 
led to instabilities during the period of shock reflection. These instabilities occurred 
because the shock waves are smeared over several mesh spaces, and the flow properties 
are mixed within the smeared shock. The instabilities can be removed by using the speed 
of sound behind the reflected shock immediately after the shock impingement in the sta- 
bility criterion instead of aoo. 


RESULTS 


The numerical method presented in this report is used to investigate the diffraction 
of a normal shock wave by a sphere and the subsequent transient phenomena in an inviscid, 
perfect-gas flow field. The transient behavior in the region ahead of the sphere is studied 
for ranges of the incident -shock Mach number M s and the specific-heat ratio y. Time 
histories are presented for the shock detachment distance and the stagnation-point pres- 
sure and tangential velocity gradient in terms of the nondimensional time parameter 


r 



(17) 


where t is the time after shock impingement, c is the speed of the reflected shock in 
the corresponding one -dimensional shock reflection problem, and 6g is the shock 
detachment distance for steady flow. 
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The flow field under consideration has three regions: the undisturbed region ahead 
of the incident shock, the region of uniform flow behind the incident shock, and the shock 
layer between the reflected shock wave and the body. 


The local Mach number in the region of uniform flow behind the incident shock is 
related to y and M s by the equation 

• i \ 2 


'i - 


Moo 2 = 


Me 


y(y - 1) 



Thus, Moo has a limiting value as M s — °° for each value of y. This limiting value 
is finite for y > 1. 

Lighthill (ref. 16) generalizes the one -dimensional shock reflection problem in 
order to treat shock reflection at a plane surface which is slightly inclined to the incident 
shock wave. To the first order in the angle of inclination 6 the expressions for the 
speed of the reflected shock wave and the thermodynamic properties behind the reflected 
shock are the same as those for the one -dimensional shock reflection problem. The 
expression for the velocity component tangent to the plane surface behind the reflected 
shock is 


v = V 


OO 


3y - 1 

y + 1 



3 - y -J-V 

- 1 Ms 2/ 


(18) 


This expression may be used to evaluate the tangential velocity component at the point 
where a normal shock wave is reflected by a sphere for the period immediately after 
shock impingement. 


Detailed Time -History Calculations 

The diffraction of the incident shock wave by the sphere and the subsequent motion 
of the reflected shock wave for y = 1.4 and M s = 5.55 (M^ = 1.70) are shown in 
figure 3. The shock waves shown in the figures are positioned at the midpoints of the 
rapid rises in the pressure profiles across the shocks. The computation is initiated at 
r = -0.21 before shock impingement has occurred. The shock wave moves downstream, 
strikes the sphere, and is reflected. Initially, the reflection is regular. When the angle 
of inclination of the incident shock to the surface has increased, Mach reflection occurs. 
Examples of regular and Mach reflection are shown in figure 3(a) for r = 0.088 and 
t = 0.37, respectively. 
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This change from regular to Mach reflection is similar to the change which occurs 
in the inclined-plane shock reflection problem. It is shown in reference 17 that the angle 
of inclination of the shock wave to the plane surface has an upper bound above which 
regular reflection is not possible, and a lower bound below which Mach reflection cannot 
occur. These bounds are functions of y and M s , and they do not coincide; there is an 
overlap region where both types of reflection are mathematically possible. 

The present method does not determine the angle of inclination at which Mach 
reflection is initiated in the shock wave-sphere interaction problem. Bryson and Gross 
(ref. 6) show that the length of the Mach shock is a very small fraction of the radius of 
the sphere for some time after Mach reflection begins. The Mach shock is defined 
clearly by the present method only after its length is greater than the radial mesh 
spacing Ar. 

The computational method must be modified in order to make calculations in the 
vicinity of the Mach shock, because the alternate use of equations (12) and (15) does not 
yield stable solutions at the Mach shock. For this reason, equations (12) are used exclu- 
sively at the points in (r, 0, t) space where Mach reflection can occur. When the Mach 
shock has moved out of the computational region, the alternating procedure is used 
throughout the flow field. 

In figure 4 the shock and sonic-line locations for t = 5.28 are compared with the 
experimental results of Ladenburg, Winckler, and Van Voorhis (ref. 18) for the steady 
flow of air over a sphere at Moo = 1.70. The numerical results for the shock-wave loca- 
tion for r = 5.28 and the experimental results are in agreement near the axis where the 
numerical solution has reached the steady state. The numerical solution is still slightly 
transient at points away from the axis and in the vicinity of the shock wave. The numeri- 
cal solution for the sonic line has the same general character as the experimental results. 
A firm conclusion cannot be drawn because of the scatter in the experimental points. The 
numerical solution does not extend to the shock wave because of the smearing of the flow 
properties in the vicinity of the shock. 

The transient stagnation streamline distributions of several of the flow properties 
for y = 1.4 and M s = 5.55 are shown in figures 5, 6, 7, and 8. These distributions are 
obtained by alternately using equations (12) and (15). Those obtained by exclusive use of 
equations (12) are compared with the results of the present method for r = 5.28 when 
the flow along the stagnation streamline has become steady. 

In figure 5 the transient pressure distributions are shown. The pressure drops 
rapidly after the initial overpressure so that the pressure profile within the layer 
approaches its steady configuration before the shock detachment distance attains its 
steady value. The pressure profile at r = 0.12 shows the maximum value of the 
stagnation-point pressure that is encountered during the computation. It is seen in the 


16 



figure that the maximum overpressure compares closely with the pressure behind the 
reflected shock in the one -dimensional shock reflection problem. The time r = 0 cor- 
responds to the time at which the incident shock should reach the body according to an 
exact one -dimensional shock computation. The numerically determined overpressure 
occurs a little later because the incident and reflected shock waves are smeared over 
several mesh spaces, and it takes several time steps for the reflection to occur. 



r/r b 

Figure 4.- Steady shock-wave and sonic-line locations for y = 1.4 and 
= 5.55 (M„ = 1.70). 

It is seen in figure 6 that the variation of the density profile with time is similar to 
that of the pressure profile. The maximum value of the numerically computed stagnation- 
point density occurs at the same time that the maximum overpressure occurs. The pres- 
ent method computes pressure more accurately than density. The numerically computed 
steady stagnation-point values for pressure and density are within 1 percent and 5 percent 
of their respective exact values. 

The distribution of the tangential velocity gradient for several times is shown in 
figure 7. It is seen that the numerical solution for the stagnation-point gradient attains 
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1.4 


1.3 



r/r b 


Figure 7.- Transient stagnation streamline tangential-velocity-gradient distributions for y = 1.4 and 

= 5-55 (M„ = 1.70). 
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its maximum value after the maximum overpressure occurs. This maximum value is 
close to the value given by equation (18). The tangential velocity gradients along the axis 
drop rapidly from the peak values and then approach the steady values slowly. 

In figure 8 it is seen that the local Mach numbers at points along the stagnation 
streamline within the shock layer adjust to their steady values even more rapidly than 
the other flow properties do. This phenomenon was observed for all the cases treated. 

The steady distributions which are obtained by using equations (12) exclusively are 
seen to be less accurate than those of the present method. The stagnation-point thermo- 
dynamic properties are more in error, and the profiles are smeared more at the shock 
wave. It should be noted that the mesh spacings which are used for the two computations 
are identical; only the methods of computation differ. 

In figure 9 the numerical results for the steady pressure distribution along the sur- 
face of the sphere for the case y = 1.4 and M s = 30 (M^ = 1.88) are compared with the 
experimental results of Kendall (ref. 19) for a sphere in air with M*, = 1.82. Kendall 
reports that his pressure measurements are accurate to within 1 percent and that viscous 
effects do not influence his results. Also, in figure 9 the numerical surface Mach number 
distribution is compared with the one given in reference 19. The results of the present 
method are seen to be in good agreement with the results of Kendall. 

Shock-Detachment -Di stance Historie s 

Numerical results for the transient behavior of the shock detachment distance for 
several cases are presented in figure 10. The nondimensional quantity 6/63 is plotted 
against r, where the values for steady shock detachment distance 6g are determined 
from the results of Van Dyke and Gordon (ref. 20) and Lomax and Inouye (ref. 21). The 
cases for y - 1.1 and y = 1.4 span a wide range of the incident shock Mach number 
M s . As Moo approaches unity the shock layer becomes large and difficult to compute 
when the steady state is approached. For this reason, the case y = 1.4 and M s = 2.51 
(M^ = 1.20) is terminated before steady flow is established. Since M*, is small for 
all values of M s when y is large, only one case is presented for y = 5/3. 

It can be seen from figure 10 that the shock-detachment-distance histories are not 
very sensitive to mesh spacing size. In this figure results are presented for two computa- 
tions for each of the cases y = 1.1 and M s = °° (M^ = 4.26) and y = 5/3 and M s = °° 
(Moo = 1.34). Finer mesh spacings were used for the computations indicated by the solid 
symbols. 

All the data points in figure 10 lie in a fairly narrow band. This result indicates 
that the nondimensional time r is an effective scaling parameter. In general, the data 
points are arranged within the band according to the value of Moo. 
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Several cases are treated for y = 1.4 to determine the asymptotic behavior of the 
shock-detachment-distance histories as M s is increased. It is seen in the figure that 
the histories approach the limiting case for fairly low values of M s ; the data points for 
M s = 30, 5.55, and 3.66 lie very near the same curve in the figure. 

An approximate analytical solution for the variation of the shock detachment dis- 
tance with time is obtained in the appendix for cases with values of y near unity and 
large values of M s . This solution for y = 1.1 and M s = °° (M^ = 4.26) is shown in 
figure 10. The constant a is evaluated by using the exact steady values of 10.00 and 
10.49 for the density ratios p s g and Pfc, respectively; the numerically obtained value 

T ? 9 

of 0.615 for vjj g, the steady nondimensional tangential velocity gradient at the body; and 

r 

the approximate value of 1 for v s , the nondimensional tangential velocity gradient at the 
shock. It is seen in the figure that the analytical solution and the numerical solution are 
in good agreement. Equation (A6) yields a value of 0.067 for 6g, as opposed to the 
numerically obtained value of 0.080. 

Stagnation-Point-Pressure Histories 

The numerical results for the stagnation-point-pressure histories are presented in 
figure 11. The nondimensional pressure P s ta^/Poo is plotted against r. Initially, an 

overpressure occurs at the stagnation point. Then the pressure drops rapidly and 
approaches its steady value monotonically. The results of Butler (ref. 9) show that the 
transient stagnation-point pressure for a flat-nosed cylinder undershoots the steady value 
when y is low and Ms is high. This phenomenon is not observed in the present results 
for the sphere. The maximum overpressures which are observed are in fair agreement 
with the pressures behind the reflected shock in the one -dimensional reflected shock 
problem, as mentioned earlier. The numerical results for the steady values are in good 
agreement with the exact values. 

A comparison of the results in figures 10 and 11 shows that the stagnation-point 
pressure approaches its steady value much more rapidly than the shock detachment dis- 
tance does. The pressures are within about 5 percent of the final values at r = 1.5. At 
this time the shock detachment distances have reached only about 70 percent of their final 
values. The shock detachment distances approach the final values at t = 5 or 6. The 
asymptotic behavior of the stagnation-point densities and temperatures is similar to that 
of the pressures. 

The initial overpressures are very sensitive to the stability criterion and the size 
of the mesh spacings. If the stability criterion is exceeded, the computations within the 
shock layer immediately after shock impingement are erratic. If the mesh spacing is 
crude, the peak values are not obtained although the solution at later times is not affected 
appreciably. This sensitivity to mesh spacing is more pronounced for low values of y. 
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(a) y = 1.4. 

Figure 11.- Stagnation-point pressure histories. 



{b} y = 1.1 and y = 5/3. 
Figure 11.- Concluded. 
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Stagnation-Point Tangential-Velocity-Gradient Histories 


The numerical results for the stagnation-point tangential-velocity-gradient histories 


are presented in figure 12. The nondimensional quantity 


—M 

Vqc \d0/ st ag 


is plotted against 


r. It is seen in the figure that the initial velocity gradient is larger than the steady value. 
After shock impingement the stagnation-point velocity gradient approaches the steady 
value monotonically. In figure 13 the stagnation-point tangential velocity gradient is 
plotted against the shock detachment distance. The velocity gradient is nondimensional- 
ized with the value which is obtained by differentiating equation (18) with respect to 8, 
and the shock detachment is nondimensionalized with its steady value. For the cases 
shown the variation of the velocity gradient with the shock detachment distance becomes 
more nearly linear as Moo is increased, and the steady gradients are approximately 
60 percent of the values given by equation (18). These observations hold for all the cases 
which were treated. 


Results for two computations with different mesh spacings are shown in figures 12 
and 13. Finer mesh spacings were used for the computations indicated by the solid 
symbols. It is seen that the stagnation-point velocity-gradient histories are fairly insen- 
sitive to mesh spacing. 


CONCLUDING REMARKS 

Numerical results are presented for the unsteady flow in the shock layer of a sphere 
which has been washed by a normal shock wave. The numerical values for the overpres- 
sure in the shock layer immediately after shock impingement are in good agreement with 
the exact calculations for the pressure behind the reflected shock in the corresponding 
one -dimensional shock-reflection problem. The results show that the growth of the 
shock-layer thickness is much slower than the adjustment of the stagnation-point thermo- 
dynamic properties to their steady values. This means that the stagnation-point pressure 
is not a good indicator of steady flow. The stagnation-point tangential velocity gradient 
also approaches its steady value faster than the shock-layer thickness for all but the 
largest free-stream Mach number. 

It is shown that the histories of the nondimensional shock detachment distance 6 / 6 g 
(where 6 s is the steady shock detachment distance) can be scaled effectively with the 
nondimensional time t = ct/ 6 g (where c is the initial reflected-shock speed and t is 
the time after shock impingement). The ratio 6/63 approaches 1 for values of r of 
about 5 or 6 for a wide range of perfect-gas specific heat ratio y and incident-shock- 
wave Mach number M s . Therefore, an estimate of the time necessary to establish 
steady flow can be obtained if estimates of 63 and c are available. 
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An approximate analytical solution is obtained for the shock-detachment-distance 
histories for cases with large values of M s and values of y near unity. 
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APPENDIX 

APPROXIMATE EXPRESSION FOR TIME HISTORY 
OF SHOCK DETACHMENT DISTANCE 

A simple analytical expression can be obtained for the time variation of shock 
detachment distance for cases with values of y near unity and large values of M s . 
Under these conditions M„ will be large, and the shock layer will be close to the body 
and can be approximated by a spherical cap in the stagnation region. In addition, the 
initial overshoots in the stagnation region are small for these conditions. 

When the reflected shock wave is approximated by a spherical surface, the tangen- 
tial velocity gradient immediately behind the shock is approximated by the free-stream 
gradient dv/d0 = Voo cos 9. This gradient along the axis is simply Voo. It should be 
noted that this approximation is good for the period immediately following shock 
impingement since the gradient given in equation (18) approaches Voo for small y 
and large M s . 

The tangential velocity component in the vicinity of the axis is assumed to be of the 
form v = v’0, where the prime denotes differentiation with respect to 9. The thermody- 
namic properties, u, and v' are assumed to be functions of r and t near the axis. 

The one-strip method of integral relations of Belotserkovskii (ref. 22) is used to 
integrate equation (5), the continuity equation for flow along the axis, across the shock 
layer. The general form for the integrated equation is 

/ \ C r=r b +6 C r=r h +6 

dr + PgUglrb + 6) + 2\ pv'dr + \ pu dr = 0 

\ / J r=r - b J r=rb 

The quantities pr, pv’, and pu in this equation are approximated with linear profiles 
across the shock layer. The integrated equation is written as 

+ y (*>b " Ps)~ + Ps*s ( r b + 1 6 ) + 6 (/>b v b + Ps v s) = 0 < A1 ) 

The shock relation 




(A2) 


32 





APPENDIX 


is used to evaluate the term p s u s . The tangential velocity gradient at the body is 
assumed to vary linearly with 6/6s. This assumption is made because M m is large 
when M s is large and y is near unity, and the variation of the velocity gradient with 
the shock detachment distance becomes more nearly linear as Moo is increased, as 
shown in figure 13. The expression for v^ is written as 


v b - v 's - ( v 's - v b,S l)^ <A3) 

The quantities in equations (Al), (A2), and (A3) are nondimensionalized as follows (the 
barred quantities are nondimensional) : 




V 

* or 


p = JL 

Poo 


6=A 

r b 


t = 


*bp 


(A4) 


where u^ = -V,*,. The equation is written as 


+ Ife + p s .M + |(p s -i)a^-i 

dt V / dt 2\ b s ) dt 2\ s ) dt 


dt 


-2 


Pb + Ps, 


3 ) V s ‘ f ) 5 _p b( v s - v b, S )|^- 0 


(A5) 


This equation can be used to obtain an expression for the steady shock detachment 
ratio 6g: 

6 S =— — 1 (A6) 

P b,S V b,S + p s,S v s " 2 

The terms dp^y/dt and dp s jdi in equation (A5) are neglected because there is 
little decompression in the stagnation region when y is near unity and Ms is large. 
For example, the initial density ratio in the shock layer immediately after impingement 
is on the order of 11.00 when y = 1.1 and M s = 00 (M^ = 4.26). The density ratios at 
the stagnation point and the shock are 10.49 and 10.00, respectively, for the steady case 

j c 

y = 1.1 and M s = 00 . The term which is proportional to 6 is more than an order of 

dt 
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magnitude less than the term p b + P g - 2^^ when y is near unity and M s is 

large. Therefore, it is neglected. 

The initial speed of the reflected shock is obtained by setting 6 equal to 0 in 
equation (A5): 



Since 


Pb + P £ 


changes very little throughout the period of unsteady flow, c is sub- 


stituted for this quantity in equation (A5). The quantities ^p^ + P^jv ' s - J- and 
Pb(v- - gj j 6g are assumed to be constant for the same reason. 

When these approximations have been made, equation (A5) can be written in the 

form 


_d_/_6_ 

dT^6 s 



“&- 1 - (1 -* 




where t 


ct 

— and the quantity a, which is constant, is given by 
6 


2p b,S v b,S " ( p b,S ‘ p s,s) v s " | 
5 b,S 7 i>,S + ? s,S ? s - I 


The solution to equation (A7) is 


6 _ 1 - e‘ m 

6 S 1 - (1 - d)e~ aT 


(A 7) 


(A8) 


(A9) 
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